Duality and the vibrational modes of a Cooper-pair Wigner crystal 
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When quantum fluctuations in the phase of the superconducting order parameter destroy the off- 
diagonal long range order, duality arguments predict the formation of a Cooper pair crystal. This 
effect is thought to be responsible for the static checkerboard patterns observed recently in various 
underdoped cuprate superconductors by means of scanning tunneling spectroscopy. Breaking of 
the translational symmetry in such a Cooper pair Wigner crystal may, under certain conditions, 
lead to the emergence of low lying transverse vibrational modes which could then contribute to 
thermodynamic and transport properties at low temperatures. We investigate these vibrational 
modes using a continuum version of the standard vortex-boson duality, calculate the speed of sound 
in the Cooper pair Wigner crystal and deduce the associated specific heat and thermal conductivity. 
We then suggest that these modes could be responsible for the mysterious bosonic contribution 
to the thermal conductivity recently observed in strongly underdoped ultraclean single crystals of 
YBa2Cu306+a; tuned across the superconductor-insulator transition. 



I. INTRODUCTION 

Ground states of superfluids and superconductors are 
characterized by a sharply defined macroscopic phase and 
an uncertain total number of particles. The mathemat- 
ical expression of this phenomenon is the number-phase 
duality: the phase tp of the superconducting (or super- 
fluid) order parameter is an operator that is quantum me- 
chanically conjugate to the particle number operator N , 
leading to the uncertainty relation AiV • A(/3 > 1. A sim- 
ilar duality exists locally. Consider, for simplicity, a lat- 
tice model of a superconductor (or a superfluid) with the 
number and phase operators on site i denoted by fii and 
(fi respectively. These operators satisfy [fii,(pj\ = iSij 
which yields the local version of the above uncertainty 



Aui ■ Aipj > Si. 



(1) 



The uncertainty relation ^ has an interesting implica- 
tion for a state of matter obtained by phase disordering 
a superconductor or a superfluid. (By phase disorder- 
ing we mean destruction of the off-diagonal long range 
order by phase fluctuations in such a way that the am- 
plitude of the order parameter remains finite and large.) 
In such a phase-disordered state the local phase is maxi- 
mally uncertain, which, according to (Q, allows the parti- 
cle number to be locally sharp. Cooper pairs (or bosons) 
may thus minimize their interaction energy by forming a 
crystal. 

A precise mathematical description of this phe- 
nomenon is given in terms of vortex-boson dualityS*^ 
which maps the system of Cooper pairs (bosons) near a 
superfluid-Mott insulator transition onto a fictitious dual 
superconductor in applied external magnetic field. The 
Cooper pair (boson) crystal emerges as the Abrikosov 
vortex lattice of the dual superconductor. 

Over the past several years it has become increasingly 
clear that the theoretical ideas summarized above may 



have found a concrete physical realization in the under- 
doped regime of high-Tc cuprate superconductors. Ac- 
cording to one school of thought the ubiquitous pseu- 
dogap phcnomenoni^ can be thought of as a manifesta- 
tion of phase-disordered superconductivity in a d-wave 
channel iSiLSi2iiS This view is supported by a number of 
experimentsiiiiSiiiii which, in essence, probe the phase 
fluctuations directly. Given this evidence and the dual 
relationship between the number and the phase one is 
compelled to ask whether the expected pair crystalliza- 
tion can be observed in underdoped cuprates. The an- 
swer to this question appears to be a tentative "yes" and 
we shall elaborate this point below. 

The experimental technique of choice to probe for 
Cooper pair crystallization is scanning tunneling mi- 
croscopy (STM). Recent advances in STM provided a 
number of atomic resolution measurements of the local 
density of states (LDOS) in a host of cuprate materials. 
These reveal an intricate interplay between quasiparti- 
cles, vortices, order and disorder under the umbrella of 
high Tc superconductivityiSiiSiiLiSii2iSa 

The STM scans have been performed on large areas of 
freshly cleaved samples of Bi2Sr2CaCu208+5 (BiSCCO) 
and Ca2-2;Naa;Cu02Cl2 (Na-CCOC) cuprate supercon- 
ductors. Real space maps of LDOS show clear modula- 
tions which can be systematically studied in the Fourier 
space. Results of the various measurements can be sum- 
marized as follows. Superconducting samples exhibit en- 
ergy dispersing features at energies that are small com- 
pared to the maximum superconducting gap. When the 
samples are underdoped toward the pseudogap phase new 
non- dispersive features appear. These features have peri- 
odicity close to four lattice constants but are not always 
commensurate with the underlying ionic lattice. At low 
energies they coexist with the dispersing featuresi^ and 
are present up to energies of the order of the gap maxi- 
mum, where the dispersing features no longer appear. 

The appearance of dispersive features in the supercon- 
ducting phase can be understood in the context of im- 



purity scattering of the quasiparticlesiSiiS^ The peak-hke 
structures in the Fourier transformed STM maps rep- 
resent the momentum transfer in these scattering pro- 
cesses. The latter are sensitive to the quasiparticle dis- 
persion and the coherence factor^ and thus aUow map- 
ping of the underlying Fermi surface, the gap function, 
and reflect on the nature of the ordered state. The dis- 
persive features also appear in the pseudogap phase,-^- 
indicating that the nodal quasiparticles survive the tran- 
sition. Furthermore, the similarity between the observed 
peaks in the superconducting and pseudogap phases sug- 
gests that the anomalous (off-diagonal) order parame- 
ter is responsible for the pseudogap phenomenoniS^ The 
existence of a residual superconducting order parameter 
in the insulating phase is consistent with the picture of 
phase disordered superconductivity. 

The non-dispersive features that are seen in un- 
derdoped samples indicate the formation of additional 
charge order close to and in the pseudogap phase. A cru- 
cial characteristic of these modulations is that they leave 
the low energy physics intact. This is manifested by the 
coexistence of the static charge modulations with the low 
energy dispersing features and by the V-shaped gap in 
the density of states. It is therefore difficult to imagine 
that an ordinary charge density wave in the particle-hole 
channel is responsible for the non-dispersive modulations. 
A conventional charge density wave would gap out the 
nodal quasiparticles and would therefore change the low 
energy properties of the system, in contrast to various 
experimental results. ^^-^A 

An interpretation of the static modulations in terms 
of the Cooper pair Wigner crystal (PWC) was first sug- 
gested by Chen et al^^ who also proposed a test designed 
to distinguish between PWC and ordinary charge density 
wave. A detailed theory of the PWC in the context of the 
phase fluctuation scenario was developed by Tesanovic,^ 
who applied the duality map to the case of a d-wave su- 
perconductor and demonstrated that scrambling of the 
order parameter phase by vortex-antivortex fluctuations 
indeed leads to Wigner crystallization of Cooper pairs. It 
was further argued that the observed patterns are con- 
sistent with a detailed account of the phase fluctuating 
(i-wave order parameter on the lattice with the relevant 
amount of charge carriers in the Cu-0 planesi^ Using 
different variations on this theme, AndersonSS and Ba- 
lents et alm^ arrived at similar conclusions. 

In this article we propose and investigate an additional 
manifestation of the existence of a PWC. Cooper pair 
crystallization should be accompanied by emergence of 
transverse vibrational modes, absent in the superfluid. 
These can be though of as transverse "magnetophonons" 
of the dual Abrikosov lattice. Both the superfluid and the 
PWC phase support longitudinal modes (phase mode or 
second sound in the former, longitudinal phonon in the 
latter). For charged systems these are gapped due to 
the long range nature of the Coulomb interaction and 
thus cannot be excited at low temperatures. Transverse 
modes, on the other hand, correspond to shear deforma- 



tions of the PWC and thus have the character of gapless 
acoustic modes. This is strictly true if one can ignore 
the coupling of the PWC to the underlying ionic lattice. 
If such coupling is present then the new modes acquire 
a gap proportional to the strength of this coupling. At 
temperatures smaller than the gap scale such vibrational 
modes would not contribute to thermodynamics or trans- 
port. We shall assume in the following that this coupling 
to the lattice is negligible although it is not a priori clear 
why this should be so. The key hint comes from the fact 
that at least in some cases the PWC has been reported to 
be incommensurate with the underlying ionic lattice^iiS 
This points to very weak coupling, for if the coupling to 
the lattice were strong PWC would always be commen- 
surate. We shall return to this important point in Sec. 
VIII. 

Using the formalism of the vortex-boson duality we cal- 
culate the normal modes of the pair Wigner crystal and 
note that such modes are accessible through the measure- 
ment of heat transport. Indeed, Taillefer et al,'^-^ have 
reported an appearance of a new bosonic mode, with 
T^ contribution to the thermal conductivity in YBCO. 
The onset coincides with the transition to the pseudogap 
state at the critical doping Xc- Our calculations below 
show that this additional thermal conductivity may be a 
result of the PWC lattice vibrations. 

This article is organized as follows. In the next section 
we briefly review a simplified version of the duality trans- 
formation that maps a phase-fluctuating superconductor 
onto the "dual" superconductor in an applied magnetic 
field. In section III we show that vortices in the dual 
model correspond to Cooper pairs and derive the inter- 
action potential between them. In section IV we consider 
an Abrikosov vortex lattice in the dual superconductor 
and calculate its normal modes. In section V we evaluate 
the inertial mass of Cooper pairs (the mass that deter- 
mines their kinetic energy) which is needed to complete 
the description of the normal modes. In the final sections 
we use the foregoing results to estimate the contribution 
of the PWC vibrational mode to the specific heat and the 
thermal conductivity, summarize our results and discuss 
their implications for the physics of cuprates. 



II. DUALITY AND WIGNER 
CRYSTALLIZATION OF COOPER PAIRS 

The route to a detailed understanding of the charge 
modulations in a pair Wigner crystal goes through a du- 
ality transformation. The need for the dual description 
has to do with the fact that Cooper pairs are highly 
non-local objects, most naturally described in momen- 
tum space when they form a coherent superconducting 
state. The duality transformation trades a description of 
a superconductor with fluctuating phase for a model in 
which the vortex field plays the primary role. The full 
theory, relevant to cuprates, must take into account the 
dynamics of vortices at the lattice scale in the presence of 



both Cooper pairs and nodal quasiparticles. Its detailed 
imple mentat ion is quite complicated and was given in 
Refs. J27l28| . Much of the complication arises precisely 
because of the need to describe phenomena on relatively 
short length scales, of the order of the period of the PWC, 
which is typically close to 4 lattice spacings. Since we 
shall be concerned primarily with the low-energy, long- 
wavelength properties of the PWC, we will contend our- 
selves with a simpler continuum version of the duality 
transformation^^ to which we input the correct lattice 
structure by hand. As usual, the long-wavelength physics 
will depend only on the symmetry and dimensionality of 
the problem. 

We thus consider an effective theory for supercon- 
ductivity in a single Cu-0 plane defined by the par- 
tition function Z = J Vl-^i ,'^*]exp{- jj^^ dr J d^xC/h) 
where ^I^ is a complex scalar order parameter and the 
Lagrangian density is 



c^Ik. 
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f 



um^). 



(2) 



The Greek index /i = 0,1,2 labels the temporal and 
spatial components of (2-t-l) dimensional vectors, d^ = 
{dcT,dx,dy) and c is the speed of light in vacuum. The 
parameters K^ = {Kq , Ki , Ki ) are related to the com- 
pressibility and phase stiffness, respectively. They could 
be calculated, in principle, from an underlying micro- 
scopic theory via the standard gradient expansion. In 
the present work we shall determine these parameters 
by matching to the relevant experimental quantities: Kq 
will be expressed through the Thomas-Fermi screening 
length and Ki through the London penetration depth. 

The electromagnetic vector potential A is explicitly 
displayed in order to track the charge content of vari- 
ous fields. In a more general model one could allow A to 
fluctuate in space and time and thus implement the long 
range Coulomb interaction; however we shall not pursue 
this here. C/ is a potential function that sets the value 
of the order parameter ^ in the superconducting state in 
the absence of fluctuations. 

Let us now proceed by disordering the phase of the 
field * = |*|e*^. We fix its amplitude |*| = *o at the 
minimum of the potential [/(l^tp) and retain the phase 
as a dynamical variable, 



C 
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K,, dj 



27r 



-A. 



(3) 



2^0 = CdT where 



Cd 



(5) 



is the speed of the phase mode in the super- 
conductor. The partition function becomes Z = 
/ V[^, **] exp(- / d^xC) with 



/: - 2^^ ( ^^^ 






(6) 



and K = VK^. 

Next we consider fluctuations in the phase. These 
can be decomposed into smooth fluctuations and singular 
(vortex) fluctuations, 



d X d^9s = 



(7) 
(8) 



The smooth fluctuations are curl-free and the singular 
fluctuations are related to the density of vortices through 
their curl in the temporal direction. 



{^x^e)o = 27TJ2q^S{r-r,iT)). 



(9) 



Here r^ (r) are the locations of vortex centers and qi are 
the respective vortex charges (-1-1 for vortices and —1 for 
antivortices) . 

In order to shift our point of view from the condensate 
to the vortices we decouple the quadratic term by an aux- 
iliary field, W^, using the familiar Hubbard-Stratonovich 
transformation. The resulting Lagrangian is 



1 



2n 



^ ^ 2K^^ + *m^m(9m^« - ^^m) + ^w^id^e,) (10) 

We may replace the third term by —idf^Wfj^Og (and a van- 
ishing surface term) and integrate over the smooth fluc- 
tuations 9s- This results in the constraint d^Wf^ = 0. In 
order to enforce this constraint we replace the field W by 
a (2-1-1) dimensional curl, W^ — {d x Ad)fi, and rewrite 
Eq. ^ as 



C^—{dxA,y^ + iidxAd)^-id,e, 



27r 

$0 



A,). (11) 



Integrating by parts in the last term and introducing vor- 
tex 3-current 



where X^ = {hKfj,/c)'^l and <l>o = hc/2e is the super- 
conducting flux quantum. To explicitly display the rela- 
tivistic invariance of the above theory it is expedient to 
rescale 



V J^o 



(4) 



and measure time in natural units related to length as 



jv = 7^-9 X d9y 

2,71 



(12) 



we obtain 

^ = W^ic)xAd)l + 2TnAd . j, - ^{d x A^) ■ A, (13) 



2K 



2111 



$n 



where the dot represents the scalar product in 2-1-1 di- 
mensions. Vortex current jy is minimally coupled to the 



gauge field Ad- This coupling implies that vortices act as 
electric charges for the dual gauge field. 

In order to complete our description of the system in 
terms of vortex degrees of freedoms we introduce a vortex 
field X that is related to the vortex current through 



ji;^^[x*d^x^{d^xnx] 



(14) 



We elaborate the theory by adding a kinetic and potential 
energy terms for the field x and write^^SiS 



1 



1 



C = \(d^^2mA^,)x\' + V{\x\) + —{dxA. 



2' 

2TTi 

---A-{dxAd) + CEM[A]. 



2K 



(15) 



where we have added, for the sake of completeness, the 
electromagnetic Maxwell term for A. The inclusion of 
the potential V restores the physics of short-range inter- 
actions between the vortex cores that has been neglected 
when we made the London approximation in Eq. P|l. 

We have now reached our goal of describing a phase 
fluctuating superconductor in terms of a vortex field x 
coupled to a dual gauge field Ad- If we ignore the cou- 
pling to electromagnetism, contained in the second line, 
the dual Lagrangian formally describes a fictitious su- 
perconductor coupled to a fluctuating gauge field. As 
we shall see shortly the temporal component of the dual 
magnetic field Bd = dx Ad is related to the charge density 
and, thus, the dual superconductor is generally subjected 
to nonzero magnetic field. 

At the mean field level this theory exhibits two dis- 
tinct phases separated by a second order transition. In 
a dual "normal" phase, (x) = 0, vortex fluctuations arc 
bounded and form finite loops in space-time. This is a 
phase coherent superconductor in the direct picture. In 
the dual "superconducting" phase the vortices condense, 
(x) 7^ 0, and vortex fluctuations proliferate. This is the 
pseudogap phase in the direct picture where phase fluc- 
tuations destroy the long range superconducting order. 

Beyond the mean field a third possibility exists in 
which the dual superconductor has a phase incoherent 
condensate. In this case (|x|) > but the dual phase 
coherence is destroyed by quantum melting of the vortex 
lattice. This phase is akin to the "vortex liquid" phases 
known from the studies of vortex matter i^ A special case 
is a situation when most of the dual vortices remain in 
a crystal but some small fraction melts and delocalizes 
through the sample. This is known to happen in a ther- 
mally fluctuating lattice superconductor— when the mag- 
netic flux per unit cell of the lattice is close to a major 
fraction, e.g. (i + -), with q ^ p. In the direct pic- 
ture this situation corresponds to the supersolid phase 
where the PWC coexists with superconductivity. Such 
coexistence seems to occur in samples of Na-CCOC—iS^ 
Whether such dual supersolid phase occurs in the class 
of quantum models considered in this paper is a very in- 
teresting open question which we shall not attempt to 
answer here. 



III. INTERACTION BETWEEN COOPER 
PAIRS 

The dual model describes a fictitious superconductor. 
We shall assume that this is a type-II superconductor 
and therefore vortices of the dual field can be present. 
Vortices in the dual model represent Cooper pairs in the 
direct picture. In order to see this we note that the dual 
flux is quantized in units of 1: compare the term |(9^ — 
27ri^^)x| in the dual model to |(9^-27ri$(7^Ap)*| in the 
direct model. Next let us obtain the relation between 
the real electromagnetic charge and the dual magnetic 
field. The coupling 2TTi^Q^A ■ {d x Ad) implies that the 
electromagnetic current is given by jfj_/hc = 27r$Q ^(i9 x 
Ad)fj., or 



J- = 2e{d X Ad)^. 



(16) 



Specifically, the flux in the temporal direction is related 
to the charge density, p = 2e{d x Ad)a. This means that 
a vortex in x carries the electric charge 



Q, = I p{v)dPr = 2e j {d x Ad)o(fr = 2e, 



(17) 



where the integration is over a surface containing the 
vortex. 

With nonzero density of Cooper pairs the dual su- 
perconductor is subject to dual magnetic field, Bda = 
{d X ^(i)o- If we assume that we are in the dual typc-II 
limit, a dual Abrikosov vortex lattice is formed. From 
now on we shall drop the adjective 'dual' whenever there 
is no potential for confusion. In the vortex lattice the 
magnetic field exhibits a periodic modulation with one 
flux quantum per unit cell. In the direct picture this cor- 
responds to a charge density modulation with charge 2e 
per unit cell. This is the pair Wigner crystal. Figure Q] 
illustrates the spatial variation of the order parameter x 
and the field Ed in such a dual Abrikosov lattice. 

We are interested in the interaction between vortices 
which we regard as point particles located at the phase 
singularities associated with each vortex. We consider a 
static configuration of vortices located at points r^ and 
neglect fluctuations in both x and Ad- This corresponds 
to a dual mean field approximation. We emphasize that 
this is a highly nontrivial mean field theory since it de- 
scribes the original superconductor in the presence of 
strong quantum fiuctuations. From Eq. ((T5)l the energy 
of a collection of static vortices becomes 



he 
~2 



-Emf ^ ^r I d'^r 



(V - 2^zAd)xl' + Ko\y X Adf] , 

(18) 

where, for simplicity of notation, we regard A^ as a 3- 
dimensional vector with z component always zero. We 
have also passed back to conventional units reversing the 
scaling introduced in connection with Eq.Q. 

The interaction energy is easiest to evaluate in the dual 
London approximation, where we assume a constant am- 
plitude of the order parameter |x(r)| = XOj but arbitrary 
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FIG. 1: (Color online) A schematic plot of the spatial vari- 
ation of the order parameter x ^nd the field Bd in the dual 
Abrikosov vortex lattice, ^d and A^ denote the dual coher- 
ence length and dual penetration depth, respectively. Dashed 
line represents the magnetic field 6d(r) associated with a sin- 
gle isolated vortex; Bd in a dense lattice is a superposition of 
such field profiles centered at individual vortices. 



variation in its phase a. As in the case of ordinary super- 
conductors the London approximation is vahd when the 
intervortex distance a 3> ^d the dual coherence length; 
i.e. the vortex cores do not overlap, ^d is difficult to es- 
timate reliably since it depends on V in Eq. H15() which 
has been essentially added by hand to account for the 
short-distance behavior of (real) vortices. We shall as- 
sume that £^d is of the order of ionic lattice spacing deep 
in the dual superconducting phase. Close to the transi- 
tion it diverges (as ~ Xo in the mean field theory) and 
thus we expect the London approximation to break down 
in this limit. 

To find the minimum of the energy we regard EyiF as 
a functional of A^ and a and vary it with respect to the 
gauge field 



5E^ 
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= 0, 



with result 



27rxo(Va-27rAd 



/fn- V X Bd. 



Now we can rewrite the energy as 
he 
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where 



X2 



4^'X^^O 



(19) 



(20) 



(21) 



(22) 



Arf is the dual penetration depth which in the context of 
the Abrikosov lattice characterizes the magnetic size of a 
vortex. In the direct picture it represents the size of the 
charge cloud associated with a single Cooper pair. 

If we apply the curl operation to Eq. H2Q|I we obtain 
the dual London equation 



-A^V' 



Brf + Bn 



zn(r) 



(23) 



where 



n{Y) =^(5(r-rj) 



(24) 



is the vortex density. The delta functions appear due 
to multiple valuedness of the phase in the presence of 
vortices; V x Va(r) — 2-kz^^6{y — Yi). It is clear that 
Bd = zBd- The solution of the London equation ^Z',\\ can 
be written as 

Bd{Y)^ I d?r'G{Y^Y')n{Y'), (25) 

where G(r) is a Green's function subject to 

{1-XIv^)G{y)^5{y). (26) 

The above equation has a simple solution in the Fourier 
space, 



G(k) 



1 



A^P- 



(27) 



Using a vector identity (V x B^)^ = B^ • (V x V x 
Bd) -I- V • (Bd X V X Bd) and discarding the vanishing 
surface term we rewrite the vortex energy (|21f) as 



Emf — 



he 

2i^ 



/dVBd-hA^V^Bd + Bd]. 



(28) 



With the help of Eqs. (|23|) and (|25|) . this can be recast 
as a density-density interaction 



Emf = \J d\ J d^r'n{Y)V{Y - r')n(r'), 



(29) 



where we have introduced the intervortex potential 
V{y) — {hc/Ko)G{Y). In view of Eq. if^ one can also 
write this as interaction between point-like particles lo- 
cated at Yi, 






(30) 



For future reference it is useful to give the intervor- 
tex potential in terms of physically accessible quantities. 
Specifically, we can trade compressibility Ko/hc for the 
Thomas- Fermi screening length Xtf- In conventional 
terms the latter is given b}*^ 



1 



X2 



= Aire - 



1 driQ 



d (9/z 



(31) 



where uq is the two dimensional density of particles and 
we have added the factor of 1/d to account for the layer 
thickness. We thus have hc/K^ = Ane'^X'^p/d. This al- 
lows us to write an explicit expression for the intervortex 
potential in the Fourier space as 



V{k) 



k^ + A-' 



(32) 



with the effective charge 



where Ri 



■tv? ^^j ) ^ij 



Ui — u," and 



47re^ /A 
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(33) 



Two remarks are in order. First, despite its appear- 
ance V{r) is not electrostatic in origin; e^ appears simply 
because we chose to express compressibility in terms of 
the Thomas-Fermi screening length. The interaction is 
mediated by phase fluctuations in the physical supercon- 
ductor. Second, it is useful to consider the interaction in 
real space. One obtains V{r) = {e1/2iT)ko{r/\d) where 
k[){x) is the Hankel function of 0-th order. Of interest is 
its asymptotic behavior. 



V{v) 




-r/\,^ r>Ad 



0.12, 



r< Ad 



(34) 



Thus the interaction has a finite range, of the order of the 
dual penetration depth. Approaching the transition, A^ 
diverges. The long distance behavior of V indicates that 
in this limit the interaction becomes long ranged but, at 
the same time, its strength vanishes since e^ — > 0. 



IV. NORMAL MODES IN A COOPER PAIR 
WIGNER CRYSTAL 

The problem of normal modes in a lattice with long 
range Coulomb interaction has been studied previously 
in the context of electronic Wigner crystalsi^ In two di- 
mensions, vibrations of a square lattice with centrally 
symmetric interactions between electrons contain modes 
with imaginary frequencies indicating that the lattice 
is, in general, unstable. In the dual picture this corre- 
sponds to the well known fact that while the triangular 
Abrikosov vortex lattice is stable, the square lattice rep- 
resents an energy maximum and is, therefore, unstablei^i 
It is also known that the square vortex lattice can be sta- 
bilized if the interactions posses four-fold anisotropyi^ 
Such anisotropics can arise from the d-wave symmetry 
of the order parameter— or from band structure effects 
associated with the underlying ionic lattice4£ In the fol- 
lowing we shall make an assumption that similar terms 
with four-fold anisotropy exist in our dual superconduc- 
tor. These terms will stabilize the square vortex lattice 
but will not affect our discussion of the long wavelength 
vibrational modes. 

We employ the standard formalism for lattice vibra- 
tions: we assume that vortices arc located at points 
Ti = 11.1 + Ui where R^ denotes vectors of a Bravais lattice 
and Ui are small displacements. The interaction potential 
(|30|l is expanded to second order in the displacements. 
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is the dynamical matrix. Aside from a constant Vq = 
Y^^, V{R.ij), the elastic energy can be written as 



AE = 2el{h~l2) 



with 






a/3 
ij J 






(37) 

(38) 
(39) 



Note that the intervortex potential V is defined every- 
where in space (not only on the lattice sites) and there- 
fore k is not restricted to the first Brillouin zone. In 
order to keep all momentum integration variables within 
the Brillouin zone we replace k — > q + Q , and 



f d^k-^Y] f d^q 
J Q Jbz 



(40) 



where Q is a reciprocal lattice vector (e*^'^ = 1 for any 
lattice vector R). 

We now analyze the terms /i and I2 by Fourier trans- 
forming the variables Uj, 



u,; = a 



BZ 



d^k 



e'l'^'Uk. 



(41) 



The momentum integral is over the first Brillouin zone 
and a is the vortex lattice constant which we include in 
order for Uk to maintain the dimension of length. Com- 
bining equations H36I41(I we find 
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where (dk) = (i^fc/(27r)^ and the integration extends over 
the first Brillouin zone. Performing the real space sums 
leads to 



AE = 2el fidk)uZu''_^Zgf3(k), 



(43) 



with 



z^pik) = J2 

Q 



{ka+Qa){kp + Qp) QaQp 



{k+Qr + xf 



Q' + K- 
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The normal modes of the system are related to the 
eigenvalues ^^(k) of the matrix Zap(}i)- Explicitly 
ruyUj'^ = 2e^Zi(k) where m„ is the dual vortex mass. We 
thus need to evaluate the reciprocal lattice sums indi- 
cated in Eq. (|44|l . These sums are slowly convergent and 
great care must be taken in evaluating them; specifically 
we need to employ the Ewald summation technique. The 
following derivation is similar to that given by Fetter— 
which was done in the context of a triangular Abrikosov 
vortex lattice. 

In order to proceed analytically we split Zap(\i.) — 
4'/3(k) + 41(k)with 



^S(k) 



KryK 



qK/3 



k^ + A- 



(45) 




(2'] 

and Z^0 (k) containing the sum over all Q 7^ 0. Evalua- 



-'ap 

tion of the latter is greatly simplified if we assume that 
A^ is much smaller than the smallest reciprocal lattice 



vector, Qo — "iir/a. Our estimate of A^ below will justify 
this assumption. Thus, 
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(46) 



We note that all the dependence on A^ is contained in 
Z^I (k) which is, in addition, independent of the lattice 

structure. .Z^o(k), on the other hand, depends on the 

lattice structure but is independent of A^. 

(2) 
The detailed calculation of Z^^ (k) is given in the ap- 
pendix. The result can be summarized as 

^a/3(k) = -P^ + ^aHSgf^k' - 2kgkp) (47) 



where for the square lattice d — 0.066. Eq. 147() is valid 
for long wavelengths, k <^ 27r/a. 

The normal modes are readily found through the eigen- 
values of the matrix Zc^(k), 



c^i(k) = 


Vgk, 




(48) 


c^2(k) = 


"'^''pa^k^ + Xf 


- 1, 


(49) 



where 



^d- 



(50) 



The first, acoustic mode, is transverse and its sound ve- 
locity is Vs- The second mode is longitudinal. At long 
wavelengths, A; <C Ad, the longitudinal mode is also acous- 
tic with velocity Vs (Xd/a) \jlfQ. As the intervortex inter- 
action becomes long-ranged (Ad — > 00) the latter speed 
diverges and the mode becomes gapped, which is ex- 
pected on general grounds. With the estimated A^ of 
about 10 vortex lattice constants the longitudinal mode is 



FIG. 2: A schematic plot of the acoustic transverse mode 
(dashed) and the longitudinal mode (solid) . The functions are 
plotted with Ad = 10a where a is the vortex lattice constant 
and uijri) is presented in units of Vs. 



unimportant (as we shall see it is the inverse of the sound 
velocity that enters the expressions for the specific heat 
and thermal conductivity). Also, had we retained the 
long-range Coulomb interaction mediated by the electro- 
magnetic gauge field fluctuations, the longitudinal mode 
would be gapped for any value of Ad. A schematic plot 
of the two modes is given in Fig. 13 



V. DUAL VORTEX MASS 

In order to evaluate the sound velocity of the transverse 
acoustic mode it is essential to estimate the mass of the 
dual vortex. Naively this mass is simply the Cooper pair 
mass, Inie- However, this is not the parameter 7n„ that 
enters the sound velocity expression. The difference be- 
tween the two originates from our scheme of calculating 
the normal modes. We have modeled the PWC as a sys- 
tem of point masses and springs. This analogy is useful 
but not quite suitable for the case of Cooper pairs that 
are delocalized over many lattice sites. The finite size of 
the Cooper pairs causes significant overlap between their 
wave functions and this leads to a situation that is quite 
different from the familiar case of vibrations of point-like 
ions in a solid. 

To make these considerations more concrete let us now 
estimate the effective Cooper pair size given by the dual 
penetration depth A^. According to the STM experi- 
ments the charge modulation amounts \,o v ~ \ — 5% of 
the total charge densityn'*^ In a dual vortex lattice this 
corresponds to the rms amplitude variation of the dual 
magnetic field v w ((S^ — -Bd)^)^/^/i3d, where Bd is the 
average field and the angular brackets denote the spatial 
average. The above rms variation depends on the ratio of 



Xd to intervortex spacing a42i Namely, a^/A^ sa z/, which 
leads to the estimate 



Ad » (5 - 10)a. 



(51) 



Near half filling we can estimate a ~ v^ag « 5A , where 
the ionic lattice constant oq = 3.8A in YBCO. Thus, 
even deep in the dual superconducting phase, the charge 
and mass of the Cooper pair are distributed over many 
lattice sites. 

In the present case it is intuitively clear that only the 
fraction of the Cooper pair mass associated with the pair 
density wave should contribute to the vibrational degrees 
of freedom. In particular as A^ — > oo the charge distribu- 
tion becomes homogeneous, the system is superfluid and 
cannot support any transverse modes. We thus seek a 
mass parameter associated with the kinetic energy of a 
moving dual vortex. The appropriate parameter is the 
inertial mass of a dual vortex. In order to determine 
the latter we note that, if we discard the coupling to 
electromagnetic field, the vortex dynamics are given by 
a relativistic theory Eq. (|15|l . Correspondingly, the rest 
energy of a dual vortex is given by 



Ey, 



m^c^ 



(52) 



Here, rriy is the inertial mass of the dual vortex that we 
seek and c^ is the dual speed of light, i.e., the phase 
velocity of the dual gauge field defined in Eq. Q. We 
may therefore estimate the energy of a vortex line fi'rest in 
the standard waji^ and deduce its inertial mass through 
Eq. ^. 

We consider a single vortex located at the origin. Its 
energy can be expressed by combining Eq. (|28|l with the 
dual London equation if^ and n(r) = S{r) as 
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bd{0). 



(53) 



Henceforth we shall denote magnetic field associated with 
a single isolated vortex by bd (r) . The energy of a vortex 
depends linearly on the magnetic field at its center. For 
a circularly symmetric vortex the London equation can 
be solved by the Hankel function 



bdir) 
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2nXl 



In 
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(54) 



where the last equality holds for r <$:, Xd- The divergence 
as r — > is unphysical and occurs due to our neglect 
of the order parameter amplitude variation in the vortex 
core. A reliable estimate is obtained by evaluating 



bd{0) - bdiU) 



1 



2nXl 



In(Kd) 



(55) 



where Kd — Xd/£,d is the dual Ginzburg-Landau parame- 
ter. Kd is assumed larger than unity (dual type-II regime) 
but since it appears inside the logarithm its exact value 
is unimportant. 



The dual vortex mass is thus given by 



incXjKi 



In(Kd) 



4e2 



MKd). (56) 



We have used the London penetration depth. 



A" 



Aire Us 



(57) 



with the superfluid density n^ = 2^Q/d to eliminate 
Ki — h'^Q/2meC. In the above estimate one should take 
the mean field value of A, i.e. the value it would have in 
the absence of phase fiuctuations. In YBCO, we thus take 
A ~ 1000 A, the value at optimum doping. Taking Kd — 
10 and d = 12A gives m„ ~ 1.1 x 10"-'^(A/Ad)^(2TOe), 
where rrie is the electron mass. 

A more instructive way of expressing m^ is to estimate 
the mean field superfluid density close to half filling by 
Us ~ l/(2agd) and obtain 



2mp 



In(Krf) 
27r 



(58) 



As expected, when Cooper pairs are localized and Xd ap- 
proaches the lattice constant oq the inertial mass of the 
dual vortex approaches that of the Cooper pair. When, 
on the other hand, Xd ^ aQ, the Cooper pair is delocal- 
ized over many lattice spacings and the dual vortex mass 
becomes small. 



VI. INTERLAYER TUNNELING AND DUAL 
MONOPOLES 

A linear dispersion such as the one we found for the 
transverse mode, combined with Bose-Einstein distribu- 
tion for phonons, leads to specific heat C^ ^ T'^ where 
d is the dimensionality of the system. The thermal 
conductivity, within the simple Boltzmann approach, is 
K — -^CvVgi where £ denotes the phonon mean free path. 
Assuming the latter is T- independent, as is the case for 
phonons scattered by the sample boundaries, we have 
K ^ T'^. We thus arrive at a conclusion that, in order to 
agree with experiment, the PWC phonons must propa- 
gate in all 3 dimensions. 

Our theory thus far focused on purely two dimensional 
physics of the Cu-0 layers. However, it is clear that vi- 
brations of PWC in the adjacent layers will be coupled. 
There are two main sources of this coupling: (i) the pair 
tunneling between the layers represented by the inter- 
layer Josephson term and (ii) the Coulomb interaction 
between the induced charge modulations. Inclusion of 
the interlayer coupling will cause the PWC phonons to 
propagate between the layers. To complete our calcula- 
tion we must determine the associated sound velocity in 
the z direction. 

Ideally, we would like to extend our duality map to 
a system of weakly coupled layers and repeat the calcu- 
lation of the normal modes. Unfortunately, there is no 



straightforward generalization of the vortex-boson dual- 
ity to a system in 3-1-1 dimensions; 2-1-1 dimensional sys- 
tems are special in this respect. The reason for this can 
be seen most clearly by returning back to Eq. (|1U|I . In 
(3-|-l)D it is not possible to enforce the zero-divergence 
constraint on VF^ by expressing it as a curl of a gauge 
field; the curl operation is only meaningful in 3 dimen- 
sions. On a more fundamental level in (3+l)D vortices 
cannot be regarded as point particles, rather they should 
be thought of as strings. Thus, rather than a dual super- 
conductor, in (3-1-1 )D the duality map produces a string 
theory. On physical grounds we still expect the PWC to 
form in a (3-|-l)D phase disordered superconductor but 
the underlying mathematical structure of the dual theory 
appears to be more complicated and beyond the scope of 
this paper 4Si^ 

We thus continue describing the physics of the individ- 
ual layers by the (2-|-l)D duality and consider the effect of 
weak interlayer coupling on the resulting quantum state. 
We focus first on the Josephson couphng, which medi- 
ates tunneling of Cooper pairs between layers. Formally 
we imagine generalizing our starting Lagrangian ^ into 
a Lawrence-Doniach mode)^^*^ for a multilayer system 
by attaching a layer index m to the scalar field ^ and 
adding a term 
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(59) 



to couple the layers. The cross terms describe tunneling 
of physical Cooper pairs between the layers. The ampli- 
tude for this process, Jo, is related to the c-axis London 
penetration depth Ac by 



Jo = 



h^c" 



d 167re2A2' 



(60) 



From the point of view of an individual layer removal 
of a Cooper pair represents a monopole event. A tunnel- 
ing event occurring at a particular instant of imaginary 
time adds or removes a Cooper pair from the plane. In 
the dual description, this event corresponds to the ap- 
pearance or disappearance of a vortex; the magnetic flux 
lines associated with it originate or terminate at the same 
point, as illustrated in Fig|3^. A point in space-time 
which represents a drain (source) of a magnetic flux is 
known as magnetic monopole (antimonopole). The vor- 
tex must reappear in the adjacent layer and this corre- 
sponds to an antimonopole. In ordinary superconductors 
vorticity is strictly conserved which can be regarded as 
a consequence of the absence of monopoles in the real 
world. In the dual superconductor vorticity is conserved 
only if we consider a system with fixed number of Cooper 
pairs. Once we introduce the Josephson tunneling vor- 
ticity is conserved only globally (i.e. the total number 
of vortices in all layers is constant) and we must permit 
monopole-antimonopole pairs to occur in the adjacent 
layers. 

To model the interlayer tunneling we should thus re- 
gard the dual gauge field Ad as compact and append to 
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FIG. 3: (Color online) a) A schematic illustration of the flux 
lines in the presence of a monopole. b) Josephson tunneling 
with lateral displacement of the dual vortex. The shaded 
regions represent charge densities associated with a Cooper 
pair in the PWC. 



the dual Lagrangian H15|l a term describing monopole- 
antimonopole events occurring in the adjacent layers. 
Since, ultimately, we only need the mean field descrip- 
tion of the intervortex interaction, we shall adopt a sim- 
pler and physically more intuitive approach. In essence, 
all we need to complete our description of PWC vibra- 
tions is an effective potential, akin to V{ri — Vj), that 
will describe the coupling between vortices in adjacent 
layers induced by the Josephson tunneling. Taking our 
clue from Eq. H3U|) we describe the system with many 
layers by a Hamiltonian H — Hq + Hj with 

^0 = ^E / / V'i(r)Vi.(r')^(r-r')Vr„(r')^™(r), 



^•' = E / / J(.r~r')[^lir)^,n+i{r') + h.c.], (61) 



The operator V'm('') creates a vortex at point r in the 
TO-th layer and J — J d^r. Hq describes interactions 
between vortices within a layer and Hj generates the 
Josephson coupling; J is the amplitude for interlayer tun- 
neling to be discussed shortly. For now we regard vortices 
as infinitely heavy (no kinetic energy in the plane). 

Our strategy will be to treat Hj as a small perturba- 
tion on the eigenstates of Hq. This is justified as long 
as J is very small. In the limit of infinite mass the un- 
perturbed eigenstates are labeled by the positions {r,i,„} 
of N individual vortices within each layer. The energy 
of this state is simply given by Eq. H30() . Clearly, there 
will be no correction to the energy to first order in Hj. 
The leading contribution appears in second order and 
describes a virtual hop of a Cooper pair from one layer 
to the next and then back. For simplicity we consider 
a case of two layers, m — 1,2, and evaluate the second 
order correction to the energy from such processes: 
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In an obvious notation \N, N) denotes a state with N 
pairs in each layer, located at {vn} and {ri2}, while 
\N + ris,N — ?]s') represents a state with one pair added 
(removed) at s in layer 1, and one pair removed (added) 
at s' in layer 2, for 77 — 1(— 1). 

The matrix element in the numerator is easily seen 
to be just |J(s — s')p but the energy denominator re- 
quires some careful thought. An obvious contribution to 
the energy difference comes from the interaction between 
Cooper pairs expressed through Eq. (|^ . For 77 = +1 
this is simply the interaction energy of a vortex added at 
s in layer 1 and antivortex added at s' in layer 2, 



AU = - J2[Vis - r,i) - V{s' - r,2)]. 



(63) 



For 77 = —1 the overall sign is reversed. However, this 
cannot be the whole story since Eq. (|63|l does not account 
for the fact that the state with TV vortices per layer cor- 
responds to the absolute minimum of the total energy of 
the unperturbed system. At the level of Eq. H3U|I this 
requirement is implemented as a constraint. This is ad- 
equate as long as N is held constant. Once we allow 
for number fluctuations, however, we must consider the 
energy cost Uq of removing (adding) a Cooper pair from 
(to) a layer. Ultimately, this cost is related to the electric 
charge neutrality: removing or adding a pair from/to a 
neutral layer costs Coulomb energy. Thus, we estimate 
Uq as the electrostatic energy of the charge distribution 
corresponding to an extra pair at s in layer 1 and a miss- 
ing pair at s' in layer 2. The missing pair is modeled 
simply as an effective positive charge distribution on the 
otherwise neutral background. We have 
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with R = (r, z) a three dimensional vector, e is the 
dielectric constant which reflects the polarizability of the 
insulating medium between the Cu-0 layers. Its value is 
around 10 in YBCO.^^ Taking 



p{R.) = pi{r- s)S{z) - pi(r - s')S{z - d). 



(65) 



where pi{r) = {2e)bd{r) is the planar charge density as- 
sociated with a single Cooper pair in PWC, we obtain 
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(66) 
with 1 = s — s'. For general values of 1 and d the above 
integral must be evaluated numerically. However, it turns 
out that we shall only require its leading behavior in the 
limit d, |1| <C Ad, i.e. the situation when the distance 
between the two charge clouds is small compared to their 
lateral size given by A^. In this limit we can expand 
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The energy correction H62() thus becomes 
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where we have suppressed various arguments in order to 
display the structure of the result. We notice that AE^^^^ 
depends on the positions {rn} and {ri2} through AC/. 
Upon integration over s, s' this will produce interactions 
between vortices in different layers, as expected. The 
origin of this interaction lies in the fact that the energy of 
the virtual intermediate state depends on the position of 
the extra Cooper pair relative to the pairs already present 
in that layer. Since the positions of the former before 
and after the tunneling event are strongly correlated this 
induces interaction between Cooper pairs in the adjacent 
layers. 

We can make the form of the interlayer interaction 
more explicit in the limit \AU\ <C Uq by expanding 
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The first term is a constant, but the second term, when 


expanded with help of Eq. (1^, contains the following 


expression, 
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(70) 
which provides an explicit interaction potential between 
the pairs located in the adjacent layers. Although we 
derived Eq. l(7n| assuming just two layers it is obvious 
that it generalizes to a multilayer system. 

In order to complete our computation of Vz we must 
determine the form of the tunneling matrix element 
J(s — s'). In the direct picture, modeled by a Lawrence- 
Doniach type Hamiltonian,*- Cooper pairs are assumed 
to tunnel 'straight up' (or down), i.e. from a point s in 
one layer to the point s in the adjacent layer. A subtlety 
in the dual picture arises from the fact that Hj describes 
tunneling of dual vortex cores and not Cooper pairs di- 
rectly. We recall that even when the position of the dual 
vortex is sharply localized at a point s, the associated 
Cooper pair charge (and number) density is delocalized 
over the length scale A^ around s with the probability 
density given by bd(j—s). Thus, we can think of a Cooper 
pair as being described by a wave function whose enve- 
lope varies as \/bd{r — s). As illustrated in Fig. |3|b a 
straight up tunneling of a Cooper pair in general may 
lead to dual vortex core tunneling with nonzero lateral 
displacement. The associated amplitude to tunnel from 
point s to s' will be given by the overlap of the two pair 
wave functions, ^ J y/bd{s — r)bd{r — s'). This last in- 
tegral is somewhat difficult to evaluate because of the 
square root, but the resulting function is clearly close to 
bd{s — s'). In the following we thus use 



J(r) « Jobd{r), 



(71) 
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which turns out to be properly normahzed since 

An important iniphcation of the above result H71|l is 
that tunneling over lateral distances larger than A^ is 
exponentially suppressed and we may indeed use the ap- 
proximation H67() for Uo when evaluating Vz from Eq. 
()70|l . In fact we will use only the first term in (|67|l which 
gives the leading contribution of the form, 
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(72) 
where we have used the relation V{y) ~ e^A^6d(r) in 
order to express everything in terms of bd- It is worth 
noting that in Eq. (|72|l the prefactor 
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of the order of the Coulomb potential at distance d, 
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Here 2i/e ~ 2eo^/A^ is the fraction of Cooper pair charge 
that is modulated. As mentioned above ;/« (1 — 5)% in 
the STM experiments. 

For the typical values of various parameters listed be- 
low Eq. 1H2) and Ad = 20A we find Ej/Ec « 30; 
the Josephson term dominates. It is also important to 
note that Ej and Ec exhibit very different scaling with 
Xd- In particular Eq. H73|) implies that Ej ^ A^ while 
Ec ^ X^ ■ Thus, even if far away from the transition the 
direct Coulomb contribution is significant, as the transi- 
tion is approached and A^ becomes large, the Josephson 
coupling between the layers always takes over. In the 
following we shall thus focus exclusively on the latter. 



characterizes the energy scale for the interlayer coupling 
while the expression in the square brackets represents a 
dimensionless function of spatial variable r whose range 
is set by A^- The interaction is repulsive. 

We could now repeat the calculation of the normal 
modes for the layered system. This calculation,*^ how- 
ever, is lengthy and does not yield any new physical in- 
sights to the problem at hand. All we need, in fact, 
is the result for the speed of sound in the z direction. 
This can be readily estimated in analogy with Eq. H50() 
which states that up to a numerical constant the speed of 
sound squared equals the energy scale of the interaction 
divided by mass. This of course is a well known result for 
ordinary systems of springs and masses; what Eq. H50|l 
confirms is that this result remains valid even in the case 
of interactions with finite range. We thus obtain 
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The above is a perturbative result and remains valid 
only as long as Jq/Uq ^ 1. With help of Eqs. (|^ and 
(|60|l this condition becomes 
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(75) 

where a = e^ /Tic ~ 1/137 is the fine structure con- 
stant. The latter appears because the Josephson cou- 
pling connects two "dual worlds" through the real world 
coupling. Since Ac is typically very large in cuprates 
- 10" - lO^A Eq. (Cni should be well satisfied except 
very close to the transition. 

We now briefly discuss the direct Coulomb interaction, 
which will turn out to be negligible in most cases. If 
we assume that only the interaction between neighbor- 
ing layers is important then standard treatment leads to 
acoustic dispersion for the transverse mode along the z 
direction with the sound velocity ^jEc/rn^ where Ec is 



VII. SPECIFIC HEAT AND THERMAL 
CONDUCTIVITY 

Having found the eigenmodes of the PWC as well as 
the effective mass of dual vortices, we may now com- 
bine these results and calculate the specific heat and the 
thermal conductivity associated with the transverse vi- 
brational modes. As mentioned above the longitudinal 
modes cither have much higher speed or are gapped and 
thus will not contribute at low temperatures. 

Let us first estimate the transverse in-plane sound ve- 
locity H50|) in terms of physically meaningful parameters. 
Combining Eqs. (|33|l and (|56() this becomes 
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We observe that the dependence on Ad has dropped out, 
except through Kd. We expect Eq. l(77| to be valid only 
deep in the PWC state; as one approaches the transition 
point the dual London approximation ceases to apply. 

Thomas-Fermi screening length is normally of the or- 
der of inverse Fermi wave vector fci?, or several A. The 
ratio of the two length scales is thus about 10~^ — 10~^. 
For Kd — 10 the square root is 0.30. Eq. (|77|l thus 
gives Vs about an order of magnitude smaller than the 
Fermi velocity in cuprates (the latter is ~ 10~^c). This 
makes sense physically since the PWC vibrations are 
purely electronic phenomenon, and, thus, on dimensional 
grounds one expects Vs and vp to be of similar order of 
magnitude. 

Similarly we can express the interplane sound velocity 
Vz associated with the Josephson coupling. Evaluating 
Eq. ^ with help of Eqs. (ESI, I®, EOI) and (gZl we 
obtain 
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Thus, Vz grows with Xd- According to the criterion 175() 
the perturbation theory breaks down when [Xd/XcOi) be- 
comes of the order of 10. A glance at Eqs. H77|l and 178|l 
confirms that this makes sense physicahy for in this hmit 
Vz approaches v^. The sound velocity becomes isotropic 
and we may no longer treat the interlayer tunneling as a 
perturbation. 

Repeating the standard calculation of the phonon spe- 
cific heat for the case of sound velocity with uniaxial 
anisotropy leads to 
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where 7 = /g dxx^ /{e' — 1) = 7r^/15. The associated 
thermal conductivity in the Cu-0 plane becomes 



with 



« = -C.vsi = /3T3, 
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Combining Eqs. (jTTj), ^^ with Eq. (|SU we obtain 

where the lower bound obtains if we take A = 1030A and 
Ac = 6.2A as found for optimally doped YBCO,"''" while 
the upper bound obtains for A = 2020A and Ac — 37. 7A, 
values relevant to the 56K "ortho-II" phase. In both 
cases we take Xtf — a-o = 3.8A , d — 12A , e = 10 and 
Krf — 10. We assume that the phonon mean free path 
t ~ 0.05mm, which is the geometric average of the width 
(^ 500^m) and the thickness (~ lO^m) of samples used 
in the experiment^ 

The experimental resulls^ for the point farthest from 
the transition gives /? « 0.5 mW/K^ cm. With the value 
of Xd = 2b — 50A, estimated from the STM data via 
Eq. H51|l . our result above is broadly consistent with the 
experimental data, although it has to be noted that it 
depends strongly on the assumed values of the input pa- 
rameters, most notably on the a5-plane and c-axis pen- 
etration depths. These are known in YBCO with high 
accuracy— but they also vary strongly with doping and it 
is not entirely clear which values one should adopt in the 
estimate of (3. We have argued above that one should take 
the underlying mean-field, noninteracting values, which 
are presumably most accurately approximated by the val- 
ues measured near the optimal doping. Within the class 
of phase fluctuation models considered here, any reduc- 
tion of superfluid density upon underdoping is attributed 
to an interaction effect beyond mean field theory. In real- 
ity, part of the change may be associated with the change 
of the underlying mean-field ground state but it is im- 
possible to make a precise statement about this. It is 



also possible that the measured values at optimal doping 
already reflect a fair amount of fluctuations and the un- 
derlying mean- field superfluid density should be higher. 
The range of values displayed in Eq. 1)82(1 is meant to be 
indicative of these various uncertainties. 



VIII. SUMMARY AND OPEN ISSUES 

We have studied the normal modes of a pair Wigner 
crystal employing a duality transformation from a phase 
fluctuating superconductor to a fictitious dual type-II su- 
perconductor in applied magnetic filed. Vortices in the 
dual superconductor represent Cooper pairs in the orig- 
inal model and the vibrational modes of the PWC can 
be calculated as magnetophonons of the dual Abrikosov 
vortex lattice. Assuming that that pinning of the dual 
vortex lattice to the underlying ionic lattice is negligible, 
as suggested by the fact that PWC is incommensurate 
in some cases, the transverse magnetophonon is acoustic 
and will thus contribute to the low energy thermal and 
transport properties of the system. For charged systems 
longitudinal modes will be gapped. Our main result is 
the estimate of the sound velocity of the transverse mode 
which determines the magnitude of its contribution to the 
specific heat and thermal conductivity. The in-plane ve- 
locity, Vs^ is found to be about an order of magnitude 
lower than the Fermi velocity vp- Together with our 
estimate of the interplane sound velocity, Vz, which we 
assumed to be determined primarily by the Josephson 
coupling between the neighboring layers, our considera- 
tions yield a T^ contribution to the thermal conductivity 
with a prefactor consistent with the recently reported 
bosonic mode in strongly underdoped single crystals of 
YBC0.-21 

An important length scale in this problem is the dual 
penetration depth, A^, which has the physical meaning 
of the size of the Cooper pair in the PWC. STM ex- 
periments indicate that A^ is much larger than the dis- 
tance between the Cooper pairs; the latter are extended 
objects with strong zero-point motion and overlapping 
wave functions. This is the main reason why duality is a 
useful concept in this problem: it provides a convenient 
tool for the description of a strongly quantum fluctuat- 
ing system of Cooper pairs in terms of dual vortices that 
can be treated as point particles. The latter are local 
objects and their physics can be accurately described in 
the mean field approximation. 

A key assumption, underlying all our preceding consid- 
erations and results, is that the pair Wigner crystal is es- 
sentially decoupled from the ionic lattice. As mentioned 
in the Introduction transverse modes of a PWC pinned to 
the ionic lattice would be gapped and thus irrelevant at 
low temperatures. The duality transformation reviewed 
in Sec. II shows that PWC indeed can exist in contin- 
uum, independently of any underlying crystalline lattice. 
There is, therefore, no logical contradiction implied by 
the above assumption. In the context of cuprates one 
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must ask to what extent does the continuum model reflect 
the physics of Cooper pairs moving in the copper-oxygen 
planes. The key issue here is whether PWC is commensu- 
rate or incommensurate with the underlying ionic lattice, 
since incommensurate PWC cannot be pinned, except by 
disorder. Pinning by disorder affects the magnetophonon 
mean free path but in general does not open a gap in the 
phonon spectrum. 

The problem of commensurability is a difficult one to 
analyze theoretically as it involves the details of PWC 
energetics, band structure and electron-ion interaction. 
Our argument in favor of the incommensurate PWC is 
therefore largely phenomenological and is based on the 
following three observations. First, the checkerboard pat- 
terns in at least some experimentsiiiS have been re- 
ported to have period clearly different from 4 ionic lattice 
constants (the values range between 4.2 — 4.7ao) imply- 
ing incommensurate PWC. The mere existence of such 
an incommensurate PWC indicates that coupling to the 
lattice must be extremely weak even in the case when 
the PWC is commensurate. Second, the PWC appears 
to exist for a relatively wide range of dopings. If only 
commensurate PWC were allowed then one would ex- 
pect dramatic changes in PWC structure upon variation 
of the pair density; in particular the unit cell size would 
vary significantly as the PWC adjusted to different dop- 
ing levels. No such dramatic variations are observed. 
It appears, instead, that PWC structure, where it ex- 
ists, is largely independent of doping. Finally, we found 
that Cooper pairs in a crystal proximate to a supercon- 
ductor have wavefunctions delocalized over many lattice 
constants. We may expect that a potential with ionic lat- 
tice periodicity should be relatively ineffective in pinning 
such delocalized objects. 

An issue which might require further consideration is 
the dependence of the sound velocity of the PWC on the 
doping level of the system. Experimentally the strength 
of the T^ mode goes to zero continuously as the tran- 
sition to the superconducting state is approached. In 
our theory the transition is marked by the divergence of 
the dual penetration depth A^- We found the in-plane 
sound velocity to be independent of A^. This is due to 
the exact cancellation between the strength of intervor- 
tex interaction e^ and the dual vortex mass my: both 
vanish as A^ — + cxd. The c-axis sound velocity is propor- 
tional to A^ which ensures that the thermal conductivity 
of the bosonic modes decreases as the transition is ap- 
proached, in agreement with the experiment. It is im- 
portant to emphasize, however, that there is no reason 
to expect that this result will remain valid very close to 
the transition. First, v^ is obtained in a perturbation 
theory which is only valid as long as Xd/^c ^ 1- Outside 
of this regime one must treat the full 3D system. Sec- 
ond, the dual London approximation also breaks down 
near the transition since the vortex cores begin to over- 
lap. To address the physics of magnetophonons near the 
transition one would have to treat vortex vibrations in 
the full Ginzburg-Landau theory Eq. (|18|) . Qualitatively 



we expect that the amplitude fluctuations will make the 
intervortex interaction stronger, leading to increase in Vs 
and thus reduction of thermal conductivity as the transi- 
tion is approached, in accord with experiment. Detailed 
calculations, however, present a daunting challenge and 
are left for future investigation. 

Another open issue is the inclusion of the detailed 
structure of the PWC. We have assumed a simple square 
Bravais lattice of Cooper pairs. ExperimentsiSiiSiSS and 
detailed theoretical considerations2i2& indicate a lattice 
with square symmetry, but with a more complicated in- 
ternal structure. Within our approach this could be mod- 
eled as a square vortex lattice with a basis. The vibra- 
tional spectrum of such a lattice will be more complicated 
but will retain the acoustic mode derived above which re- 
flects the center of mass motion of the unit cell. 

STM studies also indicate the presence of domain walls 
and other defects in the PWC. Such defects would scat- 
ter PWC phonons and cause a short mean free path £. 
Our estimate, on the other hand, suggests that in YBCO 
£ should be of the order of the sample size. This im- 
plies that the PWC in YBCO is much more homogeneous 
than that in BiSCCO and Na-CCOC. Given the extreme 
purity of the YBCO crystals used in the thermal con- 
ductivity measurement^ this is perhaps not surprising. 
Unfortunately YBCO is not amenable to high-precision 
STM studies due to its lack of a natural cleavage plane. It 
would be interesting to see if the bosonic mode could be 
observed in the thermal conductivity of Na-CCOC. Based 
on our model we would expect it to be much weaker due 
to much shorter mean free path £. 
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APPENDIX A: EVALUATION OF Z"f 



This appendix follows the calculation of Fetter f^ ad- 

(2) 

justed to the case of the square lattice. The matrix Z^n 
is a symmetric rank-2 tensor. In terms of the vector k it 
can therefore be written as 



.(2)/ 



kakfj 



Z):>{\^)^A{k^)5^p + B{k^) ^2 



(Al) 



where A and B are scalar functions of k^ . The trace of 
Z'^^ vanishes to first order in l/(A^Qo)i thus 

2A(fc2)+B(fc2) = 0. (A2) 

We evaluate the off-diagonal elements of the matrix to 
determine B{k'^). 
Let us define 



S'(k) = 27rn„ ^ 



Q#o 



\Qx + kx)\Qy + ky) QxQy 



(Q + k)2 



Q2 



(A3) 



where n„ = 1/a^ is the density of vortices. Ewald sum- 
mation technique sphts the above sum into two parts, 
S = 5'rf(k) + S'r(k), in such a way that Sd converges 
rapidly in real space while Sr converges rapidly in recip- 
rocal space. One obtains, for any Bravais latticei^ 

Sd = 2^(l-e^'^-^0^(l + '™.i?^2)e-'^"^«? 

i 3 



tjf 



2^^ y^ (Qx + K){Qy + fcy) ^-(Q+k)V47i-n^ 



27m< 






(1 



(Q + k) 



^ — fc /47ra^N 



(A4) 



where YLj = [Xj, Yj) is a lattice vector. The above sums 
can be easily evaluated numerically for arbitrary k. For 
small wave vector |fc| we may expand Sd and Sr to second 
order 



Sd 



Of 



E' 3 3 



-^^(l + 7rn,i?2)e— .H, (A5) 



-n^x'^y 



Q 






-Q^/47rn„ 
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where ^ denotes a summation that excludes the zero 
vector. In order to perform the sums for the square lattice 
we substitute 



27r 
'R = a{l,m), Q = — (l^'tn) 
a 



(A6) 



and arrive at 



l,m 



The sums over m and I are rapidly convergent and can be 
evaluated numerically. This gives the ofF-diagonal part of 

Z(2) 



Zg)(k) = -0A15^k,ky = -■da'k^k 



27r 



(A7) 



with d = 0.066. Thus, B{P) = -^a^P. 
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